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We present a novel, generally applicable Monte Carlo algorithm for the simulation of fluid systems. 
Geometric transformations are used to identify clusters of particles in such a manner that every 
cluster move is accepted, irrespective of the nature of the pair interactions. The reject ion- free and 
non-local nature of the algorithm make it particularly suitable for the efficient simulation of complex 
fluids with components of widely varying size, such as colloidal mixtures. Compared to conventional 
simulation algorithms, typical efficiency improvements amount to several orders of magnitude. 

PACS numbers: 05.10.Ln, 61.20.Ja, 64.60.Ht, 82.70.Dd 



The Monte Carlo (MC) method constitutes an impor- 
tant simulation technique in many areas of physics and 
chemistry 0, Q . One of its central aspects is the possi- 
bility to introduce non-physical dynamics, permitting the 
study of systems that evolve over otherwise prohibitively 
large time scales. A well-known example is the cluster al- 
gorithm for lattice spin models introduced by Swendsen 
and Wang (SW) 3], which suppresses dynamic slowing 
down near a critical point. Since the conception of this 
method, its generalization to off-lattice fluids of interact- 
ing particles has been an elusive goal, the main bottle- 
neck being the absence of particle-hole symmetry. Also 
away from the critical point the existence of several dif- 
ferent time and length scales constitutes a major obstacle 
in the simulation of complex fluids. This situation com- 
monly arises in multi-component systems, such as binary 
mixtures, colloidal suspensions and colloid-polymer mix- 
tures, and has essentially precluded the computational 
study of many such systems. In this Letter, we present 
a novel, rejection- free cluster Monte Carlo method of 
considerable generality that alleviates this problem. It 
greatly facilitates the canonical simulation of large classes 
of continuum systems, such as complex fluids, by gener- 
ating particle configurations according to the Boltzmann 
distribution, without suffering from severe slowing down 
in the presence of large size differences. 

The SW lattice cluster algorithm and its improvement 
by Wolff 4] are based upon the Fortuin-Kasteleyn jsj 
mapping of the Potts model onto the random-cluster 
model, which decomposes a system of spins (or Potts 
variables) into independent clusters. This is mani- 
festly different from collective update schemes in which 
more or less arbitrary groups of spins (particles) are 
flipped (moved). While such multiple-particle moves 
have yielded significant improvements in specific situa- 
tions 0, 05 Is, 9], their acceptance rate generally is an 
exponentially decreasing function of the number of parti- 
cles involved. By contrast, the SW algorithm is rejection- 
free: every completed cluster is flipped without an addi- 
tional evaluation of the resulting energy difference. Ef- 
ficient off-lattice cluster algorithms are only known for 



the Widom-Rowlinson and Stillinger-Helfand models for 
fluid mixtures 10, JJ-], in which identical particles do 
not interact at all. Furthermore, for hard-sphere fluids 
Dress and Krauth have proposed a cluster algorithm 
based upon geometric operations, that is capable of re- 
laxing size- asymmetric mixtures 13] and model glass for- 
mers 14J. However, hitherto no general off-lattice equiv- 
alent of the SW approach has been developed. In this 
paper, we demonstrate that the geometric method can be 
formulated as a Wolff single-cluster algorithm and subse- 
quently be extended to arbitrary pair potentials between 
the constituents, while maintaining its rejection- free na- 
ture. The resulting generalized geometric cluster algo- 
rithm (OCA) handles interactions in the same manner 
as the SW algorithm, and in this respect can be consid- 
ered as its counterpart for continuum systems. Lattice 
cluster methods as well as the OCA exploit invariance of 
the Hamiltonian under global symmetry operations. 

In the original OCA a molecular configuration 
is rotated around an arbitrarily chosen pivot and over- 
laid with its original (non-rotated) version. Objects that 
overlap between the original and rotated configurations 
lead to "clusters" of particles. The particles belonging 
to these clusters are exchanged between the original and 
the rotated configuration. Since the pair potential is ei- 
ther zero or infinity, each configuration without parti- 
cle overlaps has the same Boltzmann factor and hence 
the same probability. It has been suggested jf3, 
to extend this approach to other pair potentials by in- 
troducing a Metropolis-type criterion for the acceptance 
of a cluster move. However, this approach faces serious 
drawbacks, (i) It cannot be applied to soft-core poten- 
tials, since the cluster-construction process fails to gener- 
ate configurations containing interpenetrating potentials, 
(ii) The efficiency strongly deteriorates, since the algo- 
rithm is no longer rejection- free. Indeed, the stronger the 
interactions, the less relevant the (athermal) clusters be- 
come, and for many practical cases a conventional single- 
particle Metropolis algorithm will be more efficient jig . 
It is noteworthy that for lattice spin models the original 
OCA has successfully been generalized to include attrac- 
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FIG. 1: Two-dimensional illustration of the interacting geometric cluster algorithm. Light and dark colors label the particles be- 
fore and after the geometrical operation, respectively. The small disk denotes the pivot, a) Initial configuration; b) construction 
of a new cluster via point reflection of particles 1-3 with respect to the pivot; c) final configuration. 



tive short-range interactions However, this approach 
implicitly exploits particle-hole symmetry by relying on 
a mapping between sites in the original and the rotated 
lattice structure. 

In order to formulate a geometric cluster algorithm 
for interacting fluids, we first rephrase the cluster con- 
struction for the original GCA as follows. After a ran- 
dom pivot has been chosen, the first particle is picked 
at random and moved via a point reflection with re- 
spect to the pivot (this reflection replaces the rotation). 
If this leads to one or more overlaps, the correspond- 
ing particles are also moved with respect to the same 
pivot. This procedure is reiterated recursively until no 
more overlaps are present. For the next cluster, a new 
pivot is chosen. In the presence of a general pair po- 
tential F(r), we generalize this scheme as illustrated in 
Fig. n After the first particle i has been moved from 
position to its new position r-, two classes of par- 
ticles are identified: (a) particles that interact with i 
in its original position; (b) particles that interact with 
i in its new position. Every particle that belongs to 
category (a) or (b) is subsequently considered for inclu- 
sion in the cluster, i.e., for reflection with respect to the 
pivot. Particles that fall into both categories are con- 
sidered only once. While the first particle i is always 
moved, subsequent particles j are added to the cluster 
with a probability Pij = max[l — exp(— /3A^j), 0], where 
Aij = V{\t[ - r,|) - V{\r, - r,|) and P = l/k^T. Thus, 
the cluster addition probability for particle j solely de- 
pends on the energy difference corresponding to a change 
in relative position of i and j. Other energy differences 
resulting from a move of particle j are not taken into 
account, which distinguishes this method from a stan- 
dard Metropolis algorithm with multiple particle moves 
and makes it the analog of the Wolff cluster algorithm. 
Instead, the procedure is carried out iteratively. If par- 
ticle j is added to the cluster, then all its interacting 
neighbors (both in category (a) and in category (b)) that 



have not yet been added to the cluster are considered for 
inclusion as well. The cluster construction is completed 
once all interacting neighbors have been considered. 

A particle j in category (a) that is added to the cluster 
can be viewed as "moving with" particle z; a particle j in 
category (b) is then interpreted as "moving from" parti- 
cle i in its new position. However, in either case particles 
i and j maintain their original separation. The system 
evolves by virtue of the particles that are not included 
in the cluster. Again, this is analogous to a spin clus- 
ter algorithm, in which all spins within a given cluster 
maintain their relative orientation. In the limit of a pure 
hard-core repulsion, category (a) particles do not exist 
and the addition probability for category (b) particles is 
unity. Thus, the original GCA 12] indeed constitutes a 
special case of the generalized GCA. 

The ergodicity of this algorithm follows from the fact 
that there is a non-vanishing probability that a cluster 
consists of only one particle, which can be moved over an 
arbitrarily small distance, since the location of the pivot 
is chosen at random. Despite the presence of a variable 
addition probability and the existence of two categories 
of particle moves, the proof of detailed balance proceeds 
in a similar way as for the Wolff algorithm y. In the 
transition from a given configuration X (energy Ex) to 
a new configuration Y (energy £^y), an energy change is 
induced by every interacting particle that is not added to 
the cluster. The probability of such a "broken bond" k is 
given by 1— p/c, which is unity if the pair energy decreases. 
This is the analog of a pair of antiparallel spins in a lat- 
tice cluster algorithm. For an increase in pair energy A/e, 
the probability of breaking a bond is exp(— /3A/e). Ac- 
cordingly, the creation of a certain cluster corresponds to 
breaking a number of bonds, which has a probability 
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where the product runs over the set {k} of all broken 
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FIG. 2: Comparison between a conventional molecular dy- 
namics calculation (solid lines) and the geometric cluster algo- 
rithm (symbols), for a size-asymmetric binary Lennard- Jones 
mixture. Shown are, from left to right, the correlation func- 
tions for small-small, large-small, and large-large pairs. The 
agreement is clearly excellent. 



bonds, which is comprised of the subset {/} of broken 
bonds / that lead to an increase in pair energy and the 
subset {m} of broken bonds that lead to a decrease in 
pair energy. The factor C accounts for the probability 
of creating a specific arrangement of bonds inside the 
cluster. The probability for the reverse move runs over 
the same set {/c}, but all energy differences have changed 
sign (indicated by p) and the sum over A/ is replaced by 
the negative sum over the complementary set {m}, 



T(Y ^X) = C\[(l-Pk)=Ce^^ 



(2) 



The probability of picking a specific particle as the start- 
ing point for this cluster is identical in the forward and 
the reverse move. Moreover, we require the geometric op- 
eration to be self-inverse. For clusters thus constructed, 
we then indeed have succeeded in fulfilling detailed bal- 
ance while maintaining an acceptance ratio of unity: 
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Figure 121 illustrates the agreement between the general- 
ized GCA and a conventional NVT molecular dynamics 
simulation for a binary Lennard- Jones mixture contain- 
ing 800 small and 400 large particles at a total packing 
fraction r] ~ 0.213. The respective particle diameters are 
(Jii = 1.0 and <J22 = 5.0 and the interaction strengths 
equal en = 0.40 and 622 = 0.225, supplemented by the 
Lorentz-Berthelot mixing rules 2] . The particles are con- 
tained in a cubic cell with periodic boundary conditions. 
All interactions are cut off at 3<J22- Evidently, the GCA 
is capable of handling soft-core potentials. 

The true advantage of the generalized GCA transpires 
upon consideration of its efficiency. As a simple model 
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FIG. 3: Efficiency comparison between a conventional local 
update algorithm (open symbols) and the generalized geomet- 
ric cluster algorithm (closed symbols), for a size-asymmetric 
binary mixture of Yukawa particles. As opposed to the local 
algorithm, the autocorrelation time per particle (expressed in 
/is of CPU time) for the GCA depends only weakly on size 
ratio a (variations correspond to changes in the volume ratio 
of large vs small particles in the cluster), resulting in an ef- 
ficiency improvement of several orders of magnitude already 
for moderate a. 



system with intrinsically slow dynamics we again con- 
sider a binary fluid mixture of A^i small and A^2 large 
spherical particles with size ratio a = (122 /cii > 1- The 
particles are contained in a fixed volume, at equal pack- 
ing fractions ?7i = 772 = 0.1. A^2 is fixed at 150 and Ni 
increases from 1 200 to 506 250 as a is varied from 2 to 15. 
While pairs of small particles, as well as pairs involving 
a large and a small particle, act like hard spheres, the 
large particles have a Yukawa repulsion. 



U22{r) = I 



+00 r < (722 

J exp[-K(r - Cr22)]/(r/cr22) r XT22 , 



(4) 



where pj = 3.0 and the screening length = an. The 
Hamiltonian describing the system is given by the sum 
over all pair interactions. As a measure of efficiency we 
consider the integrated autocorrelation time r for the en- 
ergy 18]. For conventional MC calculations, r rapidly in- 
creases with increasing a, because large particles tend to 
get trapped by particles belonging to the smaller species 
(this situation will further deteriorate in the presence of 
an attraction between large and small particles). Indeed, 
for a > 7 it was not even feasible to accurately estimate 
r within a reasonable amount of CPU time. By contrast, 
the generalized GCA has an autocorrelation time that 
only weakly depends on the size ratio, as illustrated in 
Fig. ini At a = 7 the resulting efficiency gain already 
amounts to more than three orders of magnitude. 

To explore the performance of our algorithm near 
a critical point, we have simulated the one-component 
Lennard- Jones fluid with a potential cutoff Tc = 2.5cr at 
T* = ksT/s = 1.19 and p* = pa^ = 0.3197, very close 
to criticality . The energy autocorrelation time r ex- 
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FIG. 4: Energy autocorrelation times r vs linear system size 
for a critical Lennard- Jones fluid, in units of particle sweeps. 



hibits a power-law dependence on the linear system size L 
for local moves as well as cluster moves, see Fig.^ Since 
the density is a conserved quantity, hydrodynamic slow- 
ing down must be anticipated [2^. Owing to field mix- 
ing ^21], this will also manifest itself in the critical energy 
correlations. Thus, the acceleration ~ L^-^ achieved by 
the GCA cannot unequivocally be ascribed to the (par- 
tial) suppression of critical slowing down. This is consis- 
tent with the observation that the average relative cluster 
size grows faster than L^/^. Given the number of factors 
that determine the cluster-growth process, care must be 
taken to generalize these observations to the performance 
of the GCA at the critical point of other fluids. 

Suppression of critical slowing down is the primary 
benefit of cluster algorithms for lattice systems, making 
it a crucial requirement that the percolation threshold 
of the cluster-formation process coincides with the crit- 
ical point. The generalized GCA, on the other hand, 
addresses a much larger class of problems by accelerat- 
ing fluid simulations over a wide range of temperatures 
and packing fractions. Its essential limitation is that the 
clusters must occupy only part of the system. The av- 
erage cluster size depends not only on the interaction 
strength, but also on the total packing fraction and size 
and shape of all constituents. Although no unique perco- 
lation threshold can be defined in a continuum system of 
interacting particles, we have observed that the average 
relative cluster size increases abruptly above a certain 
packing fraction, rapidly lowering the computational ef- 
ficiency. This packing fraction depends on temperature 
and system properties, but r] ^ 0.23-0.25 represents a 
typical threshold. For increasing size ratio or degree of 
polydispersity, we expect the window of accessible pack- 
ing fractions to grow, in accordance with the increase 
of the percolation threshold as a function of polydisper- 
sity 22]. Indeed, for the binary mixtures of Fig. 01 the 
relative cluster size rapidly decreases with increasing a 
at fixed total packing fraction. 

In summary, we have introduced the first general 



rejection-free cluster algorithm for off-lattice systems. 
Its premier significance lies in a performance increase of 
many orders of magnitude for complex fluids in which 
the constituents exhibit a large size asymmetry, thus en- 
abling the simulation of mixtures that were hitherto only 
accessible via an effective one-component approach. Size 
ratios up to 100 have been reached in simulations in- 
volving several millions of particles. Our approach can 
be extended in several ways, including the treatment of 
non-spherical particles and electrostatic interactions. 
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